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A  METHOD  OF  SELF-ORGANIZING  MESH  WITH  A  COMPUTER 

Zeng  Jirong,  You  Yiren,  Shao  Yuhua  and  Liu  Tang 
Computing  Center,  Academia  Sinica 

Abstract 

In  this  paper,  we  present  a  method  of  self -organizing  mesh 
with  the  computer.  It  is  a  method  of  simple  logic  in  common  use. 
In  two  dimension  space,  this  method  produces  triangular  meshes 
on  arbitrarily  shaped  bodies  by  treating  bodies  as  collections 
of  quadrilateral  regions. 

This  method  has  been  written  in  Fortran  IV. 

The  problem  of  self-organizing  mesh  with  a  computer  is  import¬ 
ant  in  common  finite  element  software.  In  recent  years,  many 
methods  have  appeared  both  domestically  and  abroad  and  each  of 
these  methods  has  its  own  advantages  and  application  range.  The 
method  proposed  in  this  paper  is  a  method  of  simple  logic  in  com¬ 
mon  use  which  can  be  applied  to  relatively  complex  regions  with 
various  mediums  and  complex  connected  regions.  At  the  same  time, 
we  also  give  a  method  of  nodal  point  numbering.  This  nodal  point 
editing  mode  is  used  for  the  coefficient  matrix  produced  in  the 
finite  element  method  which  is  the  well  arranged  block  diagonal 
matrix,  for  the  nonzero  element  collection  as  well  as  for  band 
width  automatic  minimizing  under  certain  conditions.  After  the 
actual  application  of  electric  and  magnetic  field  computer  pro¬ 
grams,  we  obtained  better  results  than  those  obtained  abroad  with 
the  same  type  of  software  and  its  use  was  more  convenient. 

I.  Outline  of  Method 

1)  We  divided  the  region  to  be  organized  according  to  geo¬ 
metric  shape,  medium  distribution  and  requirements  into  certain 
large  "quadrilaterals"  and  set  up  a  logical  mesh  of  quadrilaterals 
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based  on  the  characteristics  of  the  boundary  lines  of  the  quadri¬ 
laterals  . 


2)  We  selected  parameters  and  and  transformed  the 

quadrilaterals  into  unit  squares  on  the  (£,»))  plane.  We 
organized  the  unit  squares  based  on  the  logical  mesh  and  obtained 
nodal  point  coordinates  (  %  *1^)  on  the  unit  squares. 

3)  We  transfored  nodal  point  coordinates  (  *1^)  on  the 

unit  squares  into  nodal  point  coordinates  {x^,y  ) . 

4)  We  put  these  quadrilaterals  together  to  carry  out  nodal 
point  numbering  and  obtained  the  nodal  point  coordinates  and  de¬ 
sired  information  on  the  entire  region. 


II.  Parametric  Transformation  and  Peripherical  Fragmenting 


It  is  assumed  that  the  parametric  equations  of  the  four  sides 
of  the  quadrilateral  are  separately 


MO,  o <# <  T„  *  — 


1,2,3, 4. 


(0 


> 


In  Fig.  1,  nu  represents  the  partition  value  required  for  side  i 

and  we  arrange  that  m2*in4  *  TaIcin9  any  one  natural  number  E, 

parameter  T.  is  divided  into  M.=*E  .  parts,  and  0=t  ^  t,  ^  t-*** 
i  i  mi  o  i  6 

S  We  calculate  the 'cumulative  chord  length  S.{t),  that  is 

Ml  1  1 

5,(0  -  2  UM*>)  “  *<(*(-.))'  +  (M*<)  -  yXii-x))1}* 

+  {(MO  —  XiOOy  +  (MO  “  M<*))'>‘/*»  '*+»  (2) 

We  introduce  parameters  ^  and  *1  and  cause 

f  *.C0  -  SM/SX  T,) ,  .  -  1,3, 

l  4,(0  -  S,(0/S,(T,)t  .'-2,4.  c 

We  let  %  and  *)  be  the  coordinates  on  the  two  orthogonal  axes  and 
thus  change  the  quadrilateral  ABCD  into  the  unit  square  A'B'C'D1 
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on  the 


(  %  ,  *1 )  plane  (see  Fig. 
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Fig.  1 


2)  . 


For  convenience,  we  eliminate  lower  symbol  i  if  we  do  not  find 
any  misreadings.  We  carry  out  partitioning  of  ArBT,  B'C',  C*B*  and 
D*  A •  based  on  the  pattern  of  change  of  the  nodal  points  on  the 
given  periphery  and  obtain  the  coordinates  of  the  peripherical 
nodal  points  on  the  (5  ,  )  plane.  The  pattern  of  nodal  point 

change  on  each  side  can  be  equidistant,  arithmetically  rise  (or 
lower)  and  geometrically  rise  (or  lower) .  The  selection  of  the 
pattern  can  be  determined  but  it  is  not  required  that  the  pattern 
of  each  side  be  the  same.  For  example,  if  the  first  side  takes 
equidistant  partitioning,  then  the  second  side  can  select  geo¬ 
metric  or  arithmetric  partitioning  and  each  side  can  even  select 
different  partitioning  modes  for  different  sections. 


III.  Generating  Logic  Mesh 


On  the  periphery,  we  can  obtain  nodal  point  coordinates 
^xi'^i^  on  the  or^9^nal  quadrilateral  periphery  from  the 
(  if  0^  coordinates  of  the  peripheral  nodal  points  based  on 
the  parametric  transformation  formula.  We  calculate  chord  length 
S j  situated  between  the  corresponding  nodal  points  on  the  two 
"vertical"  sides, 

St  —  {(*,(*,)  —  *«(»,) )J  +  (y> 0t)  —  y*(h) /  —  o,  1,2,- 


Naturally,  S0=*AB, 
ontal’ 


S_  ssCD. 
mesh  line  is  recorded  as 


The  partition  number  of  the  1  "horiz- 


n£ 


(n  «m, ,  n  *m, ).  Then 

i  IBj  J 
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We  record  the  coordinates  of  the  k  nodal  point  on  the  j£  "horizon¬ 
tal"  mesh  line  as  (k ,JL)  and^oktain  the  logical  components  of  the 
nodal  points  of  the  quadrilateral  regions.  For  example,  the 
black  dots  in  Fig.  3  give  the  logical  mesh  of  a  quadrilateral 
region.  The  outermost  layer  of  the  square  brackets  in  formula 
(4)  represents  the  taken  integer. 


H-H-r- 


i 

I 


-4-  f-  -  - 


Fig.  3 


IV.  Nodal  Point  Coordinates 


We  will  now  use  the  logical  mesh  and  nodal  point  coordinates 
(xi'^i>  on  the  plane  needed  to  produce  the  mesh  on  the 

unit  squares  of  the  (  %  ,  l) )  plane  to  realize  the  organization  of 
quadrilateral  ABCD .  For  this  reason,  we  assume  that  the  logical 
coordinates  of  nodal  points  is  and  record 


4. 

b 


(5) 
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The  corresponding  £  coordinate  is 

f.a)  -  #.«*.)>  +  +  1)  -  #,([*,  J)}(*,  -  H,D,  «  -  1,3,  (6; 

The  terms  within  the  square  brackets  represent  the  taken  integers. 
The  corresponding  nodal  point  coordinates  (  ^  (k,j£)  ,  i)(k,A))  in 
the  unit  square  A'B'C'D'  can  be  given  by  the  following  formula: 


*(*.  0 


SM*)  ♦  (Mk>)  -  tMM.  0. 

nXJ)  +  (*0)-*(0)UU 

i -(%0)  -  w(0X&(4»)*-C«.)>' 


(7) 


(*.  y> 

After  obtaining  £  (k,j2)  and  *J(k,jG),  we  use  interpolation^ to  find 
nodal  point  coordinates  (x(k,<£),  y(k,jf)),  that  is 

”  f  *(*,  o  -  oo-  o)  +  «uot,  /))fo.  o, 

‘  y(A*0  —  *($(♦»  0)0  —  *0*0)  +  y*(f0»0)n0»0,  (  ' 

In  the  formula,  jt=»l,2***,  m2-l,  k*l,2,.**,  n^. 

V.  Correction  of  the  Nodal  Point  Coordinates 

In  order  to  improve  the  shape  of  the  element  near  the  boundary 
for  the  "quadrilaterals"  composed  of  certain  boundary  curves,  we 
can  introduce  parameters  4  and^  in  the  calculations  to  correct 
%  and  and  thus  improve  the  shape  of  the  element.  Below  we 
will  consider  two  different  correction  methods. 

1)  To  correct  boundary  parameters  and  (i-1,3,  j=2,4) , 

we  use  ^  and  separately  to  replace  and  Formula 

(2)  can  be  written  as 

“  StCif-0  +  AJj(7V)»  (2') 

In  the  formula 
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A5XT,)  -  AS ,(T,)(l  +  (i  “  2)”ga(*'  ""  2'}*)’  ‘  “  1’3* 

AS'(ry)  -  as,(t,) (i  +  |^y;j’|  (3  -  0  -  2*>f)>  «  -  2>  4> 

AS4(Tf)  —  S£Ti)  -  5.(7*,-,),  «  —  1 , 2,  3,  ♦. 


The  meaning  of  4y.  and  4x.  is  similar  to  ^S..  In  the  form- 

1  1  1 

ulas,  4*  and  &  are  determined  by  the  properties  of  the  boundary 
curves.  In  >rder  to  simplify  the  calculations,  can  be  taken 
as  a  natural  number,  ^  as  a  real  number  larger  than  zero,  for 
example,  letting  *2  and  /3  **  1. 

2)  We  used  /)  ^  (£)  and  (Z)  separately  in  the  calculations 
of  formula  (7)  to  replace  *)  2  (fi)  and  04 1*)  and  obtain 

*(0  -  *<0  -  +  3  ~  0  - 1.(0) 

AJ, (I)  ' 

X  #Hgn(Ar,(i))»  »  “  2,  4,  (9) 

In  the  formula 

A*.(0  -*<(/  +  !)-*.(/-  1), 

Ay,(0  -  *0  +  1)  -  *(l  -  l), 

AS,(0 -(A,,(0*  +  AyXO*)*^ 

In  the  above  formula,  4*  and  ^  are  determined  by  the  properties  of 
the  boundary  curves.  In  order  to  make  calculations  more  conven¬ 
ient,  we  can  take  *  to  be  a  natural  number  and  ft  is  taken  as  a 
real  number  larger  than  zero.  Therefore,  4  S.^  (fi)  can  be  used 
to  substitute  in  the  following  formula,  that  is 

^(O-UxXOI  +  I  Ay,  (01. 

vi.  Generating  Triangular  Elements 


In  two  adjoining  rows  of  nodal  points,  each  row  successively 
takes  two  points  to  form  a  quadrilateral.  We  calculate  the  length 
of  the  two  diagonal  lines  and  join  the  two  vertices  of  the  short 
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line  to  form  a  triangular  element.  See  Fig.  4.  We  first  use  the 
four  points  of  A,B,C  and  D  to  form  quadrilateral  ABCD.  If 
BD  <  AC,  then  we  join  BD  to  form  <AABD  and  successively  take  E  to 
form  quadrilateral  BDCE.  If  BC  <  DE , then  we  join  BC  to  form 
^BDC,...  This  continues  until  the  nodal  points  of  the  two  rows 
are  completely  joined.  Therefore,  we  obtain  the  desired  triangu¬ 
lar  element. 

The  above  organization  method  is  not  only  applied  to  "quadri¬ 
lateral"  regions  but  is  also  applied  to  "triangles"  and  even  "two 
sided  figures"  (e.g.  Fig.  5).  In  Fig.  5,  A  and  B  coincide,  C  and 
D  coincide  and  it  is  only  necessary  that  the  nodal  point  number 
of  the  two  "vertical  sides"  be  the  same.  The  dividing  points  of 
the  upper  and  lower  "horizontal"  sides  can  be  arbitrary  which 
includes  the  use  of  zero. 


Fig.  4 


•I  V,  s/ 

\l4-V 


Fig.  5 


The  sides  which  can  be  defined  as  "vertical  sides"  or  "hor¬ 
izontal  sides"  are  relative.  Actually,  for  a  "quadrilateral," 
it  is  only  necessary  to  have  the  number  of  dividing  points  of  a 
pair  of  sides  be  equal.  As  regards  the  whole  calculating  region, 
it  is  necessary  that  the  side  of  this  "quadrilateral  "  be  a 
"vertical  side"  and  that  of  the  adjoining  "quadrilateral"  be  a 
"vertical  side."  Whichever  side  is  selected  as  the  "vertical 
side"  should  be  determined  by  the  properties  of  the  region's 
shape  and  solution.  The  terms  "quadrilateral",  "vertical  side" 
or  "horizontal  side"  etc.  used  in  this  paper  are  only  for  conven¬ 
ience  of  narration. 


Use  of  the  mesh  organized  by  this  method  has  the  following 
advantages:  (1)  The  boundary  nodal  points  accurately  fall  on  the 
boundary  line.  (2)  We  can  control  the  ’density  of  the  nodal  points 
according  to  demands.  (3)  The  mesh  lines  have  certain  similarities 
(shape  and  dimension)  to  the  boundary  lines  and  changes  are  rela¬ 
tively  smooth.  (4)  It  can  combine  artifical  organization  and 
self-organization.  (5)  It  can  be  applied  in  various  regions  and 
it  is  convenient  for  handling  various  types  of  mediums  and  com¬ 
plex  connected  regions.  (6)  It  is  easy  to  extend  to  three- 
dimensional  space  regions.  (7)  The  amount  of  operations  is 
relatively  small. 

VII.  Nodal  Point  Numbering 

The  nodal  points  obtained  from  the  above  method  can  use  two 
methods  for  nodal  point  numbering. 

1)  Stratified  Editing 

Each  "quadrilateral"  region  is  joined  together  and  is  sequen¬ 
tially  numbered  from  left  to  right  and  from  bottom  to  top.  This 
numbering  method  has  the  following  advantages:  1.  Permutation  is 
even,  logic  is  simple  and  editing  is  convenient.  2.  The  coef¬ 
ficient  matrix  is  a  block  diagonal  matrix  and  within  it  the  block 
on  the  principal  diagonal  line  is  also  a  diagonal  matrix.  The 
block  on  the  secondary  line  has  a  nonzero  element  concentration, 
there  are  no  zero  elements  between  the  nonzero  elements  and  the 
positions  of  the  nonzero  elements  can  be  conveniently  calculated. 
Therefore,  the  coefficient  matrix  is  suitable  for  a  condensed 
storage,  the  storage  amount  is  small  and  it  is  convenient  for 
access  as  well  as  finding  solutions.  3.  Under  certain  conditions, 
the  proper  determination  of  the  "vertical  direction"  can  cause  the 
bandwidth  of  the  coefficient  matrix  to  be  very  small  and  not  to 
require  optimization  processing.  The  limiting  factor  of  this 
method  is  that  it  is  required  that  the  dividing  point  number  of 
the  "vertical  side"  be  equal  on  the  entire  region. 


Fig.  6  Computer  organized  mesh  calculated  by  a  permanent 
magnetic  moment  machine's  magnetic  field. 

2)  Block  Numbering 

After  each  small  "quadrilateral"  is  organized,  we  immediately 
carry  out  numbering.  After  this  is  completed,  we  then  combine 
the  nodal  points  of  the  small  "quadrilaterals"  and  delete  the 
vertically  repeated  numbering  of  the  repeated  side  points  in  order 
to  cause  the  bandwidth  of  the  coefficient  matrix  to  reduce  and  to 
carry  out  optimization  processing.  This  numbering  system  does  not 
require  opposite  sides  with  equal  dividing  points  and  it  is  con¬ 
venient  for  a  mesh  in  the  vicinity  of  dense  single  points.  How¬ 
ever,  the  majority  of  advantages  of  stratified  editing  are 
possibly  lost. 

The  above  method  has  used  corresponding  FORTRAN  IV  edited 
software  which  can  be  used  on  a  computer  equipped  for  FORTRAN  IV 
language.  Actual  application  shows  that  very  good  results  were 
attained  for  the  organization  of  various  regions.  In  electro¬ 
magnetic  field  computations,  the  organization  results  were  better 


than  those  attained  abroad  on  the  same  type  of  software  and  its 
use  was  more  convenient.  Figures  6  and  7  are  two  examples.  The 
figures  were  drawn  by  computer  and  the  thick  lines  in  the  figures 
are  the  contour  lines  of  the  regions. 


Fjg.  7  Computer  organized  mesh  calculated  by  a  C  type  magnet's 
magnetic  field. 
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A  HYBRID  SCHEME  FOR  THE  COMPUTATION  OF  FLUID  DYNAMIC  EQUATIONS 

Chen  Guangnan 

Abstract 

On  the  basis  of  Godunov's  formulas  for  the  resolution  of  a 
discontinuity  with  weak  wave  approximation,  this  paper  presents 
a  self-adjusting  hybrid  scheme  for  computation  of  fluid  dynamic 
equations.  This  scheme  deals  effectively  with  both  shock  and 
rarefaction  waves  and  obtains  relatively  clear  oscillation-free 
monotonic  transition  of  discontinuity  while  maintaining  high  order 
of  truncation  error  in  the  smooth  part  of  the  solution.  In  com¬ 
puting  contact  discontinuity,  the  results  are  not  as  good  as 
those  of  the  linear  discontinuity  resolution  scheme.  This  paper 
gives  simple  explanations  of  the  stability  and  solution  of  main¬ 
taining  monotonicity  for  the  obtained  hybrid  schemes.  Two  models 
are  used  to  numerically  calculate  the  compressible  ideal  fluid 
and  comparisons  are  made  with  other  finite  difference  schemes. 

I. 

We  consider  the  one -dimensional  nonstationary  fluid  dynamic 
equations 

X"  W+-£-F(W)-Q.  (1) 

Of  Ott 

In  the  formula  IK  —  (»,•,  E)T,  F(W)  —  p,  p*)1  .  Equation  (1) 

can  also  be  written  in  the  equivalent  form 

±-W+  (2j 

At  _  ftr 

In  the  formula,  matrix  A  has  the  following  form 

o—i  o' 

*  "  dp/d*  —  wdp/dt  dp/d*  ,  (3) 

.  »dp/d0  p  —  Sdp/dt  mdp/da. . 

Here,  v,  u  and  E  separately  represent  the  total  energy  of  the 
specific  volume,  momentum  and  unit  mass.  Total  energy  E  is  the 


sum  of  the  internal  energy  and  kinetic  energy 

E  *  e  +■  —  h1, 

2 

and  between  internal  energy  e  and  pressure  p  and  specific  volume 
v  is  satisfied  state  equation 

p=p(v,e) 

We  know  from  the  thermodynamic  relational  formula  that  the 
Lagrangian  velocity  of  sound  is 

S-f  ~  df/ZS. 

It  is  easy  to  find  that  the  characteristic  value  of  matrix  A 
corresponds  to  0,  +c. 


As  regards  equation  (1),  if  we  use  Godunov's  formulas  for  the 
resolution  of  a  discontinuity  with  weak  wave  approximation  [1], 
then  we  have 

I  um+l  —  _ P-’k 


■f* 

^4  -  ^  -  p-u*> 


In  the  formulas 


I U7  -Uf*--'  - 


f?*i*?*i  +  ‘7-i*7-i  -  (P?+t  — 
<7+{  t 


L-  _  -  m’-t) 

[  P?  -  Prtt,‘  - 7  4. , - 

ef*i  + 

Here,  C  is  found  based  on  equations  (4)  and  (5). 

We  can  deduce  the  equations  which  draw  near  to  difference 
schemes  (6)  and  (7)  as 

|£-s-<*+->- 

j!r,'£(',  +  ,)-11  c») 

V  ““  ^  a)  *  0 

dr  dr  v 


1  A  dp  e  ft 

—  —  —  Ax  — q  — - Ax  — 


in  the  equations  2edx  2  dr.  •  We  can  see  by 

comparing  equations  (8)  and  (1)  that  each  of  the  equations  pos¬ 
sesses  the  smoothing  effects  of  terms  m  and  q  which  is  generally 
called  scheme  viscosity.  Because  of  the  existence  of  linear 
scheme  viscosity,  the  difference  schemes  only  have  first  order 
accuracy.  In  order  to  overcome  the  excessive  dissipation  caused 
by  the  scheme  viscosity  and  raise  the  accuracy  of  the  schemes,  we 
must  try  to  eliminate  the  linear  scheme  viscosity  and  cause  the 
difference  scheme  to  possess  second  order  truncation  errors. 


For  this  reason,  the  calculation  of  the  velocity  and  pressure 
of  the  "lattice  points"  in  equation  (7)  separately  add  on  cor¬ 
rection  terms  (  J0)n  and  (  4p)^.  Here,  Av  and  Ap  are 

1  3 

undertermined  parameters  and  we  use  them  to  cause  the  difference 
schemes  to  satisfy  second  order  accuracy.  We  substitute  equation 
(7)  with  correction  terms  Av  and  into  equation  (6),  carry 

out  Taylor  series  expansion  on  the  [  ( j+1/2)  ^x, n  At]  point,  use 
second  order  accuracy  and  then  have 


8>r 

a  , .  -  a» 

j_l 

—  Ari£.\ 

_  JL, 

ft  +  2 

dl3 

ft  + 

dx  ‘ 

ft* 

&.  +  JL 

a*. 

A  ,  + 

JM 

'iA,^) 

+  li 

ft  -  2 

dr3 

ft 

ax  ' 

'  2  ax/ 

ft 

d>E 

Ai  +  -5fiL  — 

.  _d_ 

('ax&-' 

)-A 

ft  2 

d/3 

Or 

flr 

\2f  ax  - 

1  dx 

+  -£-(«  •  AP)  +  —■  ip  •  AU)  +  0(  Ai3  +  Ax3)  -  0. 
Or  Or 


Using  equation  (1)  and  relational  formula 

_ &JL,  &1L.  -  A-(ei*L)  frE  -  a  f  ,  dm  ,  dp\ 

a*1  w  ft*  d*l  a,/’  v  ax  V  d7  *£)• 

after  arrangement  and  discretization,  we  can  find  the  expression 
of  the  correction  terms  to  be 


•  *r-i  /  2  ^  \ 

(AP)7  —  •  2  +  ej-i  ~  Ax)  7 


I 


When  these  correction  terms  are  added  to  equation  (7)  ,  we  then 
have 


u;  -  -  0;  -  1  (P%i  -  pr-t), 

2  Ar 

P?  “  P7‘t,)  “  PT  —  T  Ajt  _  “*-*) 


In  the  equations 


o;  - 


c7*i  +  <7-* 


Pr- 


+  *?-**•«■* 
*?4  +  <7-* 


do) 


In  this  way,  we  then  obtain  conservation  type  difference 
schemes  (6)  and  (9)  with  second  order  accuracy.  Matrix  A  and  its 
derivative  to  not  clearly  appear  in  the  scheme  and  therefore  cal¬ 
culations  are  relatively  simple. 


II. 


Because  the  second  order  accuracy  scheme  is  not  a  monotonic 
scheme,  it  easily  produces  oscillation  behind  the  wave  when  cal¬ 
culating  the  problem  of  discontinuity.  Therefore,  we  must  try 
to  add  an  artificial  viscosity  term  in  the  difference  scheme  so 
as  to  cause  the  equation  to  become  am  equation  with  viscous  dis¬ 
sipation.  When  we  add  artificial  viscosity  term  Q  into  equation 
(9)  ,  we  then  have 

o-  -  ur*>  -c?  -7(^  +  0*)  -  *'-*)» 

P7  -  pjw  -  v,  -  ±  +  p?)  4444(44  -  ■'-*)» 

In  the  equations,  Q*  “  y  —  c’-i\/(cT*i  •  *“-*)  .  Here,  b  is  a  dim¬ 
ensionless  factor  and  it  is  called  the-  viscous  coefficient. 
Difference  schemes  (6)  and  (11)  with  viscous  terms  still  possesses 
second  order  accuracy.  It  possesses  dissipation  effects  in  the 
discontinuity  area  and  therefore  the  oscillation  amplitude  of  the 
solution  is  controlled.  The  amplitude  shrinks  with  the  enlarge¬ 
ment  of  viscous  coefficient  b  yet  we  cannot  completely  eliminate 
the  oscillation  behind  the  wave. 
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Another  method  for  obtaining  dissipation  is  the  use  of  the 
self-adjusting  hybrid  scheme  [2] .  It  can  guarantee  the  mainten¬ 
ance  of  high  order  accuracy  in  the  smooth  area  of  the  solution, 
cause  the  scheme  to  obtain  a  sufficient  amount  of  dissipation  in 
the  discontinuity  area,  cause  the  scheme  to  maintain  monotonicity 
and  eliminate  oscillation  on  the  discontinuity. 

The  hybrid  scheme  is  a  scheme  formed  from  two  different 
schemes  based  on  certain  convex  combinations 

+  0  -9)L,}W:,  (12) 

In  the  formula,  is  the  first  order  accuracy  scheme,  L 2  is  the 
high  order  accuracy  scheme  and  dimensionless  quantity  9  is  called 
the  self-adjusting  switch  0  £  9  ^  1.  It  controls  the  calcula¬ 

tion  of  the  changes  of  the  area  and  plays  the  role  of  promptly 
transforming  and  calculating  the  scheme. 

When  we  take  difference  schemes  (6)  and  (7)  as  and  differ¬ 
ence  schemes  (6)  and  (9)  as  L 2,  we  can  then  establish  the  follow¬ 
ing  hybrid  schemes. 

^  (Utt’  -  U?V) 

y*7*i*?*.(tf*i  -#♦!>/(*&*  • 

+  -i  a?  *•(#.*  -  •  cr-i)  a*) 

2  0 

«;;;  -  -  pt*’’) 

+  ~  -  «r+0  -  j e-K-(uU(  -  *;-*)  ( H) 

-  E7+i  - — mw#?  -  p 

*  Ax 

+  i-  v.M.APr+xiriii  -  •  c’+i)  +  o?+l(«r*j  -  «.;♦*)] 

—  ^67K’[V7(p7*\  -  ?7-i)/(‘7+i  •  c7-()  +  C7(«r+i  -  «?-*)] 

-  -  —  -  *?+*)  (pm  -  *%*)/( 

4  A( 

+  ±**VK;Kr(-7+i  -  'T-i)(rT+i  -  •  <■?-*>.  0^> 

4  A* 
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In  the  formulas,  the  expression  of  U?"-*',  P7,L>'  is  (9)  and  the 
expression  of  C^.P?  is  dO)  . 


At  *  <■/-* 

Ax 

At  2e%±  •  c’-i 


—  !fM  ‘  c'-i  /  A/  c%\  + 

^  "  *-  +cf-iV  Ax  2  }' 

At  ci*i  + 


A/  -‘/♦i  *  / 


Ax 


-)- 


Self-adjusting  switch  dj  is.  taken  as 


•7  “  ^-4) 

<y  i4 
<2im. 


,  !  ~  rf+j' 1  ~  I  r’+j  -  »f-+  I  |  m 

i  •  *>;„{  —  v%±  I  +  |  *7^  —  r%j  1 1  *  "=*  '*’**♦  ” 

0. 


(16) 


(17) 

'?**!  +  l»M  -  *7-*  I  >  «*, 

(18) 


Key:  (1)  When;  (2)  Otherwise. 


Here,  *•  —  0.01ma*!xf+4  —  xf.j!  .  if  we  use  the  following  discretiza¬ 
tion  approximate  expressions: 

.  .  -  — 
f  P'+i  ~  Pi-i  -  -  *r+i*7-*(*?-4  -  ) »  .(19) 

l  E?+i  —  £?-*  -  —  P?(x?^  -  •%{)  +  W*?*i  -  «•;.*). 

and  substitute  them  into  hybrid  schemes  (13) -(15),  we  can  then 
obtain  the  following  hybrid  scheme: 


-(L,F*)-* 

4-  i-  {er+tK%,(pU *  -  p?+i)  -  e?/C(p?*i  -  «'?-*)} *  (-0) 

2 

—  (tj*»*)i+* 

+  i-  {dr+lxr+.(«"*{  -  *;♦*)  -  WK'(*?+i  -  do 

-  (£,£■>,♦* 

+  i  {Sf+1K7+I(E7+i  -  £/♦*)  -  O'KXE'+i  -  £,”-*)} 

z 

+  — ■  ~~~  {dt^Kf^K^CuM  *?♦()  (*’»+(  *’>+1^ 

4  Aj 

-  ff7K*K ••(«;♦*  -  •?-*)  (•)’♦*  -  <~) 
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Here,  the  L ^  operator  indicates  difference  schemes  (6)  and  (9). 

Based  on  research  by  Yanenko,  Shokin  as  well  as  Hirt  on  the 
stability  theory  for  hyperbolic  type  equations  and  difference 
schemes  [3],  we  can  deduce  that  the  first  differential  approxima¬ 
tion  of  Godunov's  schemes  (6)  and  (7)  with  first  order  accuracy 
is 

A.  W+  JL  F(W)  -  -  c  w\,  (23) 

Ar  8br  #r  1 2  \  Ax '  bx  J 


and  the  second  differential  approximation  of  second  order  accur¬ 
acy  schemes  (6)  and  (9)  is 


Equations  (23)  and  (24)  are  the  second  order  and  fourth  order 
parabolic  equations  respectively.  When  conditions 


are  established,  equations  (23)  and  (24)  are  both  satisfied.  To 
obtain  the  relationship  of  the  compatability  of  the  differential 
equations  and  the  stability  of  the  difference  equations,  we  must 
cause  difference  schemes  (6)  and  (7)  and  difference  schemes  (6) 
and  (9)  to  be  stable  and  must  satisfy  the  conditions  in  (25). 

Furthermore,  from  equation  (12)  we  obtain 

Di.ll  <  dllL.B  +  (l  — 

This  explains  that  when  the  and  schemes  are  both  stable, 
hybrid  schemes  (13)- (15)  are  also  stable.  Moreover,  in  the  same 
way  we  must  satisfy  the  conditions  in  (25). 

It  is  generally  considered  that  if  a  certain  operator  L  acts 


on  any  monotonic  mesh,  operation  LW  of  function  W  is  also  mono¬ 
tonic  and  it  is  then  said  that  finite  difference  factor  L  main¬ 
tains  monotonicity.  When  there  is  a  constant  coefficient  scalar 
equation,  hybrid  schemes  (20) -(22)  can  be  expressed  as 

wQ  —  (£.»), >4 

-  ( L* );+i  +  1  Kier«(~?+ 1  -  »U\)  ~  91  ( -  -r-i ) J ,  (26) 

In  the  formula, 

K  "  *  (*  ”  e  ^)»  “  ■“(*»♦*•  *>-*)«  ,  |A*«r(), 

—  mf-i ,  fl(x,  y)  -  !  *  — 1 ! . 

*  +  y. 

It  can  be  proved  [2]  that  if  the  conditions  in  (25)  are  estab¬ 
lished,  then  scheme  (26)  maintains  monotonicity. 


Therefore,  in  the  smooth  calculation  area  of  0=O(^3x), 
self-adjusting  hybrid  schemes  (20) -(22)  possess  second  order 
accuracy  but  when  0=0(1)  and  it  is  also  in  the  discontinuity 
area,  the  schemes  possess  the 


1  (a *Y  9 

2  At  9m 


term  of  dissipation.  Because  of  this,  when  calculating  discontin¬ 
uity,  the  solution  is  relatively  smooth  and  monotonic  and  at  the 
same  time  can  be  applied  evenly  and  pulled  wide.  In  order  to 
obtain  a  clear  and  steep  discontinuity  section,  it  is  of  signifi¬ 
cance  to  use  the  artificial  compression  method  [2] .  It  can 
cause  the  width  to  be  processed  into  a  certain  width  shock  wave 
using  the  yh  3  t  quantity  level  elongated  contact  discontinuity 
transition  area  and  cause  the  originally  discrete  shock  wave  to 
become  even  steeper.  Furthermore,  the  artificial  compression 
method  can  be  separated  from  the  main  calculation  process  and  be 
processed  independently. 


III. 


We  use  the  above  difference  scheme  to  carry  out  numerical 


calculations.  For  convenience  of  discussion,  we  consider  the 

ideal  fluid  and  at  this  time,  the  specific  formulas  of  equations 

2 

(4)  and  (5)  are  p=(Y-l)e/v,  c  =  Y*p/v.  In  these  formulas,  V  is 
the  specific  heat  ratio  and  we  take  Y=1.4.  In  order  to  test  and 
compare  difference  schemes  (6)  and  (7) ,  difference  schemes  (6)  and 
(9)  and  difference  schemes  (13) -(18)  as  well  as  the  processing  of 
the  artificial  compression  method,  we  selected  two  typical  models. 

Model  I.  The  problem  of  the  propagation  of  the  shock  waves. 

We  separated  one  shock  wave  into  two  constant  states,  v  *1,  u1*l, 
e^«3.929,  p.^«1.5716;  v2“2'  u2=0,  e2=2.858,  p2=0.5716.  Appended 
symbols  1  and  2  separately  indicate  the  states  behind  and  in  front 
of  the  wave.  The  shock  wave  used  in  this  test  starts  from  x=50  at 
the  initial  moment  and  its  propagation  velocity  is  equal  to  1. 

The  numerically  calculated  mesh  ratio  -  4-  ■  0.2.  Figures  1-4 
give  the  distributions  of  v>u,p  and  E  when  in  t*70At.  The 
straight  line  is  the  accurate  solution  and  the  curve  is  the  numer¬ 
ical  solution.  At  this  time,  the  shock  wave  appears  in  the 
x=64  area. 


<r 


i 


Fig.  3  Fig.  4 

Model  (I)  Discontinuity  occurs  in  the  x=50  area 

when  -  =  0.2,  t=70  ^  t,  t=0 . 

— . — . —  Linear  discontinuity  resolution  scheme 

-  Second  order  accuracy  scheme 

-  Self-adjusting  hybrid  scheme 

Artificial  compression  method 


Model  II.  We  consider  the  problem  of  the  diaphragm.  We 
assume  that  the  states  of  the  left  and  right  sides  of  the  dia¬ 
phragm  in  the  x*50  area  are  v^=2.245,  u^=0.698,  ej*19.796, 
p^*3.5272  and  V2»2,  Uj^O,  ©2*^*®®®'  P2=0.5714.  When  t=0,  the 
diaphragm  breaks  and  produces  a  shock  wave  towards  the  right  and 
a  rarefaction  wave  towards  the  left.  Contact  discontinuity 
appears  in  the  x*50  area.  Numerically  calculated  mesh  ratio 


t»70<4!  t. 


Figures  5-6  give  the  distributions  of  v  and  u  when 


Calculation  results  show  that:  1.  The  solution  of  linear  dis¬ 
continuity  resolution  schemes  (6)  and  (7)  with  first  order  accur¬ 
acy  are  monotonic  and  smooth.  However,  when  calculating  the 
shock  waves  and  rarefaction  waves,  the  transition  area  is  pulled 
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relatively  wide  and  generally  spans  6-8  mesh.  2.  We  can  maintain 
good  gradient  when  calculating  the  discontinuity  for  second  order 
accuracy  schemes  (6)  and  (9)  yet  intense  oscillation  can  occur 
behind  the  wave.  It  is  advantageous  for  second  order  accuracy 
schemes  (6)  and  (11)  with  viscosity  to  have  reduced  oscillation 
amplitude  behind  the  wave  and  the  suitable  selection  of  viscous 
coefficient  b  causes  the  amplitude  to  decrease.  However,  we  are 
still  unable  to  completely  eliminate  the  oscillation.  The  calcul¬ 
ation  results  of  schemes  (6)  and  (11)  are  not  drawn  in  the  figure. 

3.  The  self-adjusting  hybrid  schemes  (13)- (18)  are  very  complete 
in  eliminating  oscillation  behind  the  wave.  It  guarantees  that 
the  discontinuity  section  possesses  a  certain  gradient  and  also 
maintains,  in  principle,  understood  monotonicity.  The  physical 
picture  is  relatively  clear  and  smooth.  In  comparison  with  resol¬ 
ution  schemes  of  discontinuity  with  weak  wave  approximation,  the 
results  are  better  when  calculating  shock  waves  and  rarefaction 
waves;  when  calculating  contact  discontinuity,  the  transition 
area  is  wider  and  the  results  are  poorer. ^The  artificial  com¬ 
pression  method  can  cause  the  hybrid  scheme  to  be  better  processed 
when  calculating  discontinuity.  It  causes  even  more  abrupt 
changes  in  the  shock  wave  transition  area,  the  pictures  in  the 
contact  discontinuity  area  to  be  even  clearer  but  does  not  affect 
the  rarefaction  waves. 


Fig.  5  Fig.  6 

Model  (II)  =  0.2,  t*70  4  t. 

Resolution  scheme  for  linear  discon¬ 
tinuity 

Second  order  accuracy  scheme 
Self-adjusting  hybrid  scheme 
Artificial  compression  method 

Finally,  we  would  like  to  thank  comrade  Li  Deyuan  for  his 
valuable  views  on  this  paper. 
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